#!/usr/bin/python

from Bio import SeqIO
import subprocess
import re
from Bio.Emboss.Applications import NeedleCommandline
class compare:
	
	def __init__(self,fastaf1,fastaf2):
		
		# Load Fastas
		self.subject = SeqIO.read(fastaf1,'fasta')
		self.target =  SeqIO.read(fastaf2,'fasta')
		self.subject_tms = self.hmmtop(fastaf1)
		self.target_tms = self.hmmtop(fastaf2)
		self.scores = []
		self.alignments = {}
		
		# Make sure we have TMSs before we even continue
		if self.subject_tms is False or self.target_tms is False:
			return None
		self.cross_all()
		
	def hmmtop(self,fasta_file):
		cmd = "hmmtop -if=%s sf=FAS -is=pseudo -pi=spred" %(fasta_file)
		handle = subprocess.Popen(cmd,shell=True,stdout=subprocess.PIPE)
		results = handle.communicate()[0]
		tms = re.compile("(IN|OUT)\s+([0-9]+)\s+([^\n]+)?")
		line = tms.search(results).groups()
		if int(line[1]) is 0:
			return False
		keys = {}
		ranges = re.compile("(\d+)\s+(\d+)")
		for x in enumerate(ranges.finditer(line[2])):
			keys.setdefault(x[0]+1,[]).append(int(x[1].groups()[0]))
			keys.setdefault(x[0]+1,[]).append(int(x[1].groups()[1]))
		return keys
		
	def needle_by_tms(self,subject,target):
		needle=NeedleCommandline(gapopen=100,gapextend=50)
		needle.asequence = "asis:"+self.subject.seq[self.subject_tms[subject][0]-1:self.subject_tms[subject][1]]
		needle.bsequence = "asis:"+self.target.seq[self.target_tms[target][0]-1:self.target_tms[target][1]]
		needle.outfile="stdout"
		scores = re.compile("Length:\s+(\d+).+Score:\s+([\d.\d]+)",re.DOTALL)
		alignments = needle()[0]
		self.alignments[str(subject)+str(target)] = alignments
		(length,nwscore) = scores.search(alignments).groups()
		return float(nwscore)/int(length)
		
	def cross_all(self):
		for subject in self.subject_tms.keys():
			for target in self.target_tms.keys():
				data=(self.needle_by_tms(subject,target),subject,target)
				self.scores.append(data)
		self.scores.sort(key=lambda x:x[0], reverse=True)
		
	 